#! /usr/bin/env python

import os
import sys
import numpy
import nibabel
import mdp

ROOT_PATH="/data/FXTAS/CORRECTED"
#NIFTI_FILES_PATH="/home/xqwang/projects/FXTAS/CORRECTED/ntrimouse"
NIFTI_FILES_PATH=os.path.join(ROOT_PATH,"ntrimouse")
#CORRECTED_NIFTI_FILES_PATH="/home/xqwang/projects/FXTAS/CORRECTED/ntrimouse_corrected"
CORRECTED_NIFTI_FILES_PATH=os.path.join(ROOT_PATH,"ntrimouse_corrected")
#LANDMARK_FILES_PATH="/home/xqwang/projects/FXTAS/CORRECTED/landmarks"
LANDMARK_FILES_PATH=os.path.join(ROOT_PATH,"landmarks")
#CORRECTED_LANDMARK_FILES_PATH="/home/xqwang/projects/FXTAS/CORRECTED/landmarks_corrected"
CORRECTED_LANDMARK_FILES_PATH=os.path.join(ROOT_PATH,"landmarks_corrected")
#ROI_LANDMARK_FILES_PATH="/home/xqwang/projects/FXTAS/CORRECTED/ROI_landmarks"
ROI_LANDMARK_FILES_PATH=os.path.join(ROOT_PATH,"ROI_landmarks")
#MASK_FILES_PATH="/home/xqwang/projects/FXTAS/CORRECTED/masks"
MASK_FILES_PATH=os.path.join(ROOT_PATH,"masks")
#PLAYGROUND="/home/xqwang/projects/FXTAS/CORRECTED/playground"
PLAYGROUND=os.path.join(ROOT_PATH,"playground")

def distort(x,y,z):
    '''
    Given true World Coordination (x,y,z), return distorted (x',y',z').
    The distortion in x-direction (as a function of z) is 1.1839341-0.0004113*z^2. 
    The distortion in y-direction (as a function of z) is 1.0606435-0.0003936*z^2. 
    '''
    
    return (x*(1.1839341-0.0004113*z**2),y*(1.0606435-0.0003936*z**2),z/0.86)

def restore(dx,dy,dz):
    '''
    The inverse of distort
    '''
    return (dx/(1.1839341-0.0004113*(0.86*dz)**2),dy/(1.0606435-0.0003936*(0.86*dz)**2),0.86*dz)

def interpolate_tl(data, ii, jj, kk):
    '''
    (ii, jj, kk) are given as a location in terms of the voxel coordinate system, but may not be
    exactly on the lattice. The intensity value of of the given location is returned.

    Tri-linear interpolation is used to determine the intensity. 
    (An explanation of tri-linear interpolation can be found at  http://en.wikipedia.org/wiki/Trilinear_interpolation)
    
    '''
    i_f=numpy.floor(ii)
    i_c=i_f+1
    i_delta=ii-i_f

    j_f=numpy.floor(jj)
    j_c=j_f+1
    j_delta=jj-j_f
   
    k_f=numpy.floor(kk)
    k_c=k_f+1
    k_delta=kk-k_f

    try:
        Vi_00=data[i_f, j_f, k_f]*(1-i_delta)+data[i_c, j_f, k_f]*i_delta
    except IndexError:
        Vi_00=0

    try:
        Vi_01=data[i_f, j_f, k_c]*(1-i_delta)+data[i_c, j_f, k_c]*i_delta
    except IndexError:
        Vi_01=0

    try:
        Vi_10=data[i_f, j_c, k_f]*(1-i_delta)+data[i_c, j_c, k_f]*i_delta
    except IndexError:
        Vi_10=0

    try:
        Vi_11=data[i_f, j_c, k_c]*(1-i_delta)+data[i_c, j_c, k_c]*i_delta
    except IndexError:
        Vi_11=0

    Vj_0=Vi_00*(1-j_delta)+Vi_10*j_delta
    Vj_1=Vi_01*(1-j_delta)+Vi_11*j_delta

    #assert Vj_0*(1-k_delta)+Vj_1*k_delta >= 0

    return max(Vj_0*(1-k_delta)+Vj_1*k_delta,0)

def interpolate_nn(data, ii, jj, kk):
    '''
    (ii, jj, kk) are given as a location in terms of the voxel coordinate system, but may not be
    exactly on the lattice. The intensity value of of the given location is returned.

    Nearest neighbour interpolation is used to determine the intensity. 
    '''

    i_f=numpy.floor(ii)
    if ii-i_f < 0.5:
        i=i_f
    else:
        i=i_f+1

    j_f=numpy.floor(jj)
    if jj-j_f < 0.5:
        j=j_f
    else:
        j=j_f+1

    k_f=numpy.floor(kk)
    if kk-k_f < 0.5:
        k=k_f
    else:
        k=k_f+1

    try:
        return data[i,j,k]
    except IndexError: # out of bound
        return 0


def image_distortion_correction(img_id,domain=NIFTI_FILES_PATH,codomain=CORRECTED_NIFTI_FILES_PATH):
    '''
    Framework for distortion correction.
    '''

    domain=os.path.expanduser(domain)
    codomain=os.path.expanduser(codomain)

    if not os.path.exists(codomain):
        os.mkdir(codomain)

    if os.path.splitext(img_id)[1] == '': # img_id is a pure id, not a file name
        input_brain_image=os.path.join(domain, img_id+".nii.gz") # here I hard-coded the suffix filename, bad taste, fix it later 
        output_brain_image=os.path.join(codomain, img_id+".nii") # here I hard-coded the suffix filename, bad taste, fix it later 
    else:
        input_brain_image=os.path.join(domain, img_id) 
        if not os.path.expanduser(domain) == os.path.expanduser(codomain):
            output_brain_image=os.path.join(codomain, img_id) 
        else:
            output_brain_image=os.path.join(codomain, os.path.splitext(img_id)[0]+"_cor"+os.path.splitext(img_id)[1])
   

    # From nii_file read out necessary informations: dimensions, spacings, origin, orientation, etc.
    D_img=nibabel.load(input_brain_image)
    shape=D_img.shape
    dtype=D_img.get_data_dtype()
    D_hdr=D_img.get_header()
    D_data=D_img.get_data()
    i_spacing=D_hdr['pixdim'][1]
    j_spacing=D_hdr['pixdim'][2]
    k_spacing=D_hdr['pixdim'][3]
    x_offset=-D_hdr['qoffset_x'] # (-) changes back from the Nifti World coordinate to DICOM's scanner coordinate system
    y_offset=-D_hdr['qoffset_y'] # (-) changes back from the Nifti World coordinate to DICOM's scanner coordinate system
    z_offset=+D_hdr['qoffset_z']
    dicom_affine=D_hdr.get_qform()
    dicom_affine[0]=-dicom_affine[0]
    dicom_affine[1]=-dicom_affine[1]
    dicom_affine_inv=numpy.linalg.inv(dicom_affine)
    
    # Construct empty images of the given shape.
    data=numpy.zeros(shape,numpy.int16)

    K=data.shape[2] # number of slices      
    I=data.shape[0] # left/right 
    J=data.shape[1] # up/down 

    for k in range(K):
        #logger.debug("working on %d th slide", k)
        z=z_offset+k*k_spacing
        (xd,yd,zd)=distort(0,0,z)
        (id,jd,kd,a)=numpy.dot(dicom_affine_inv,[xd,yd,zd,1])
        if kd < 0: 
            continue
        elif kd > K:
            break

        for j in range(J):
            y=y_offset+j*j_spacing
            (xd,yd,zd)=distort(0,y,z)
            (id,jd,kd,a)=numpy.dot(dicom_affine_inv,[xd,yd,zd,1])
            if jd < 0:
                continue
            elif jd >= J:
                break

            for i in range(I):
                # convert the voxel position to the scanner position.
                # currently, just to make the new image have the same origin as the old image, and no tilting
                x=x_offset+i*i_spacing
                (xd,yd,zd)=distort(x,y,z)
                (id,jd,kd,a)=numpy.dot(dicom_affine_inv,[xd,yd,zd,1])
        
                if id < 0:
                    continue
                elif id > I:
                    break 
                else:
                    intensity=interpolate_tl(D_data,id,jd,kd) # Tri-linear interpolation
                    data[i,j,k]=numpy.int16(intensity)


    # construct nifti image
    
    # more header information need to be transfered, later ...
    hdr=D_hdr.copy()
    # new qform/sform in header, preserving the physics location of the first voxel (I need to reconsider this) and removing the tilting
    aff=numpy.array([[-i_spacing,0,0,-x_offset],[0,-j_spacing,0,-y_offset],[0,0,k_spacing,z_offset],[0,0,0,1]])
    hdr.set_qform(aff, code='scanner')
    hdr.set_sform(aff, code='scanner')

    img=nibabel.Nifti1Image(data,affine=aff,header=hdr)
     
    img.to_filename(output_brain_image) 

def mask_distortion_correction(mask_id,domain,codomain):
    '''
    Before the image distortion was discovered, all the brain segmentation works have been done. This utility
    correct the distortion of the segmentation masks.
    An old mask was traced out from ROI images and thus applys only on the ROI images. After this correction, we will generate a mask that
    applys on the full distortion-corrected image. More precisely, a new corrected mask should 
         -. have the same shape, spacing, offsets as the corrected image. 
    '''

    domain=os.path.expanduser(domain)
    codomain=os.path.expanduser(codomain)

    if not os.path.exists(codomain):
        os.mkdir(codomain)

    # From the mask_id construct the image file name (distorted image) 
    D_image_file=NIFTI_FILES_PATH +'/'+ mask_id+".nii.gz" 

    if os.path.splitext(mask_id)[1] == '': # mask_id is a pure id, not a file name
        input_mask_image=os.path.join(domain, mask_id+".img") # here I hard-coded the suffix filename, bad taste, fix it later 
    else: # mask_id is a file name
        input_mask_image=os.path.join(domain, mask_id) 
        mask_id=os.path.splitext(mask_id)[0]

    if domain != codomain:
        output_mask_image=os.path.join(codomain, os.path.splitext(mask_id)[0]+".nii") # here I hard-coded the suffix filename, bad taste, fix it later 
    else: # resolve file name conflict
        output_mask_image=os.path.join(codomain, mask_id+"_cor"+".nii")
         
    # From image file read out necessary informations: dimensions, spacings, origin, orientation, etc.
    D_img=nibabel.load(D_image_file)
    D_shape=D_img.shape
    dtype=D_img.get_data_dtype()
    D_hdr=D_img.get_header()
    #D_data=D_img.get_data()
    i_spacing=D_hdr['pixdim'][1]
    j_spacing=D_hdr['pixdim'][2]
    k_spacing=D_hdr['pixdim'][3]
    x_offset=-D_hdr['qoffset_x'] # (-) changes back from the Nifti World coordinate to DICOM's scanner coordinate system
    y_offset=-D_hdr['qoffset_y'] # (-) changes back from the Nifti World coordinate to DICOM's scanner coordinate system
    z_offset=+D_hdr['qoffset_z']

    # From image file read out the mask data.
    M_img=nibabel.load(input_mask_image)
    M_data=M_img.get_data()
    M_K=M_img.get_shape()[2] # number of slices      
    M_I=M_img.get_shape()[0] # left/right 
    M_J=M_img.get_shape()[1] # up/down 

    M_hdr=M_img.get_header()
    m_i_spacing=M_hdr['pixdim'][1]
    m_j_spacing=M_hdr['pixdim'][2]
    m_k_spacing=M_hdr['pixdim'][3]

    # While the the coordinates of landmarks and roi_landmarks were calculated, it was assumed that the resolution 0.075x0.075x0.25 for think-sliced image
    # and 0.075x0.075x0.5 for thick-sliced image, it was also assumed that the image from is not tilted wrt to Dicom reference frame. But either of this
    # assumptions turned out to be true. Now that we need to take this into consideration while correcting the distortion.
 
    # From the mask_id to obtain the landmark and ROI landmark position, and then compute the relative position of the mask in
    # the distorted image
      
    lm_text_filename=LANDMARK_FILES_PATH+'/'+mask_id+".txt"
    roi_lm_text_filename=ROI_LANDMARK_FILES_PATH+'/'+mask_id+".txt"
    with open(lm_text_filename,"rb") as lm_f:
        a,b,c=[float(x) for x in lm_f.readline().split()]
    with open(roi_lm_text_filename,"rb") as roi_lm_f:
        A,B,C=[float(x) for x in roi_lm_f.readline().split()]
    
    if D_shape[2] <=50: # thick sliced image
        rel_loc=[round((a-A)/0.075), round((b-B)/0.075), round((c-C)/0.5),1] # relative location in terms of voxles in thick sliced image 
    else:
        rel_loc=[round((a-A)/0.075), round((b-B)/0.075), round((c-C)/0.25),1] # relative location in terms of voxles in thin sliced image 
   
    # Compute the affine transformation that maps the voxel index in the mask to the *DISTORTED* physical coordinate 
    ###########
    # The affine transformation from the voxel index space of the mask to the *DISTORTED* phyiscal coordinate is given by
    # (i,j,k)^T |--> D_affine * (i+rel_loc[0], j+rel_loc[1], k+rel_loc[2], 1)^T 
    # where D_affine is the map from voxel index space tot he physical coordinate in the *DISTORTED* image (see code section below). 
    # The inverse transformation is then
    # (xd,yd,zd)^T |--> D_affine_inv*(xd,yd,zd,1)^T - rel_loc
    ###########
    D_affine=D_hdr.get_qform()
    D_affine[0]=-D_affine[0]
    D_affine[1]=-D_affine[1]
    D_affine_inv=numpy.linalg.inv(D_affine)

    # Construct empty images of that should have the same shape as the original distorted image 
    data=numpy.zeros(D_shape,numpy.uint8)

    for k in range(D_shape[2]):
        z=z_offset+k*k_spacing
        (xd,yd,zd)=distort(0,0,z)
        #logger.debug("z=%f, zd=%f", z, zd)
        #(id,jd,kd,a)=numpy.dot(mask_affine_inv,[xd,yd,zd,1])
        (id,jd,kd,a)=numpy.dot(D_affine_inv,[xd,yd,zd,1]) - rel_loc
        if kd < 0:
            #logger.debug("Not in mask area yet, kd=%d", kd)
            continue
        elif kd >= M_K:
            #logger.debug("Out of the mask area, kd=%d", kd)
            break

        for j in range(D_shape[1]):
            #import pdb; pdb.set_trace() 
            y=y_offset+j*j_spacing
            (xd,yd,zd)=distort(0,y,z)
            #(id,jd,kd,a)=numpy.dot(mask_affine_inv,[xd,yd,zd,1])
            (id,jd,kd,a)=numpy.dot(D_affine_inv,[xd,yd,zd,1]) - rel_loc
            if jd < 0:
                #logger.debug("Not in mask area yet, jd=%d", jd)
                continue
            elif jd >= M_J:
                #logger.debug("Out of the mask area, jd=%d", jd)
                break

            for i in range(D_shape[0]):
                x=x_offset+i*i_spacing
                (xd,yd,zd)=distort(x,y,z)
                (id,jd,kd,a)=numpy.dot(D_affine_inv,[xd,yd,zd,1]) - rel_loc

                if id < 0:
                    #logger.debug("Not in mask area yet, kd=%d", kd)
                    continue
                if id > M_I:
                    break 
 
                intensity=interpolate_nn(M_data,id,jd,kd) # nearest neighbour interpolation
                data[i,j,k]=numpy.uint8(intensity)


    # construct nifti image
    
    # more header information need to be transfered, later ...
    hdr=D_hdr.copy()
    # new qform/sform in header, preserving the physics location of the first voxel (I need to reconsider this) and removing the tilting
    aff=numpy.array([[-i_spacing,0,0,-x_offset],[0,-j_spacing,0,-y_offset],[0,0,k_spacing,z_offset],[0,0,0,1]])
    hdr.set_qform(aff, code='scanner')
    hdr.set_sform(aff, code='scanner')

    img=nibabel.Nifti1Image(data,affine=aff,header=hdr)
     
    img.to_filename(output_mask_image) 

def landmark_distortion_correction(img_id, domain, codomain):
    '''
    Domain is the directory that contains (distorted) landmark coordinates in .txt format.
    Codomain is the directory to store the distortion corrected landmark coordinates (in .txt format). 
    '''

    domain=os.path.expanduser(domain)
    codomain=os.path.expanduser(codomain)

    if not os.path.exists(codomain):
        os.mkdir(codomain)

    # From the img_id construct the image file name (distorted image), transformation file name, and the output image file name (corrected image)
    id=os.path.splitext(img_id)[0].split('_')[0]
    date=os.path.splitext(img_id)[0].split('_')[1]
    #D_image_file=NIFTI_FILES_PATH +'/'+ id+date[4:6]+date[0:4]+".nii.gz"
    D_image_file=NIFTI_FILES_PATH +'/'+ id + '_' + date + ".nii.gz"
    D_landmark_file_name=domain + '/' + id + '_' + date + '.txt'
    if domain != codomain:
        C_landmark_file_name= codomain + '/' + id+"_"+date+".txt"   
    else: # resolve file name conflict
        C_landmark_file_name= codomain + '/' + id+"_"+date+"_cor.txt"   
        
         
    # From nii_file read out necessary informations: dimensions, spacings, origin, orientation, etc.
    D_img=nibabel.load(D_image_file)
    shape=D_img.shape
    dtype=D_img.get_data_dtype()
    D_hdr=D_img.get_header()
    D_data=D_img.get_data()
    i_spacing=D_hdr['pixdim'][1]
    j_spacing=D_hdr['pixdim'][2]
    k_spacing=D_hdr['pixdim'][3]
    x_offset=-D_hdr['qoffset_x'] # (-) changes back from the Nifti World coordinate to DICOM's scanner coordinate system
    y_offset=-D_hdr['qoffset_y'] # (-) changes back from the Nifti World coordinate to DICOM's scanner coordinate system
    z_offset=+D_hdr['qoffset_z']
    dicom_affine=D_hdr.get_qform()
    dicom_affine[0]=-dicom_affine[0]
    dicom_affine[1]=-dicom_affine[1]
    #dicom_affine_inv=numpy.linalg.inv(dicom_affine)

    # Read out the landmark information
    with open(D_landmark_file_name,"rb") as D_lm_f:
        D_landmarks=D_lm_f.readlines() 

    with open(C_landmark_file_name,"wb") as C_lm_f:
        for lm in D_landmarks:
            lm=map(float, lm.strip().split())    
            # get the voxel coordinate
            id=lm[0]/i_spacing
            jd=lm[1]/j_spacing
            kd=lm[2]/k_spacing

            #import pdb; pdb.set_trace() 
            # get the MRI scaner coordinate
            (xd,yd,zd,a)=numpy.dot(dicom_affine,[id,jd,kd,1])

            x,y,z=restore(xd,yd,zd)
 
            C_lm_f.write("{0} {1} {2}\n".format(x,y,z))



def image_ortho_and_distortion_correction(img_file,domain=NIFTI_FILES_PATH,codomain=PLAYGROUND):
    '''
    1). Correct the distorted image.
    2). Place the corrected image in a common reference frame base on the PCA configuration of the geometric landmarks.
    #3). Correct the distorted mask, and place the corrected mask in the same reference frame of the corrected image.
    
    Input: A distorted image given by "img_id", found in the "domain" directory
    Output: 1). A corrected and orthognalized image, saved in the "codomain" directory. A corrected image means the distortion in the MRI is removed.
                A orthognalized image means 
               -). The center of the image coincide with the the center of the image (determined by the center of the landmark configuration with some adjustment), and
               -). The 3 PCA eigenvectors (computed from the corrected landmark configuration) conincide with the "World coordinate".
            #2). A corrected brain mask for the corrected image. (That is, the corrected mask should be the brain characteristic function for the corrected image.)   
    '''

    domain=os.path.expanduser(domain)
    codomain=os.path.expanduser(codomain)

    if domain==codomain:
        codomain=codomain+"/cor/"

    if not os.path.exists(codomain):
        os.mkdir(codomain)

    ## From the img_file construct the image file name (distorted image), transformation file name, and the output image file name (corrected image)
    image_ID=os.path.splitext(os.path.splitext(img_file)[0])[0]
    mouse_ID=image_ID.split('_')[0]
    date=image_ID.split('_')[1]
    D_image_file=os.path.join(domain,img_file)
    C_image_file=os.path.join(codomain,img_file)    
    C_landmark_file_name=os.path.join(CORRECTED_LANDMARK_FILES_PATH, mouse_ID + '_' + date + '.txt')
         
    # From nii_file read out necessary informations: dimensions, spacings, origin, orientation, etc.
    D_img=nibabel.load(D_image_file)
    D_shape=D_img.shape
    dtype=D_img.get_data_dtype()
    D_hdr=D_img.get_header()
    D_image_data=D_img.get_data()
    i_spacing=D_hdr['pixdim'][1]
    j_spacing=D_hdr['pixdim'][2]
    k_spacing=D_hdr['pixdim'][3]
    dicom_affine=D_hdr.get_qform()
    dicom_affine[0]=-dicom_affine[0]
    dicom_affine[1]=-dicom_affine[1]
    dicom_affine_inv=numpy.linalg.inv(dicom_affine)

    with open(C_landmark_file_name,"rb") as C_lm_f:
        C_lm_lines=C_lm_f.readlines()

    lm_conf=numpy.array([])
    for lm_line in C_lm_lines:
        lm=map(float, lm_line.strip().split())
        lm_conf=numpy.concatenate((lm_conf,lm))
    lm_conf=lm_conf.reshape((8,3))

    # constructing the landmark configuration:
    lm_pca_node=mdp.nodes.PCANode()
    lm_pca_node.train(lm_conf)
    lm_pca_node.stop_training()
    (z_sd, x_sd, y_sd)=map(numpy.sqrt, lm_pca_node.d)
    center=lm_pca_node.avg[0]+[0,0.8*y_sd,0.8*z_sd]  # adjust the pca center so that the brain fit into the image the best
    e_vect=lm_pca_node.get_projmatrix(transposed=0) 
    if e_vect[0][2] < 0:    # eighen vector 0 is the direction of brain length, should be aligned with the positive z direction
        e_vect[0]=-e_vect[0]
    if e_vect[1][0] < 0:    # eighen vector 1 is the direction of brain width, should be aligned with the positive x direction
        e_vect[1]=-e_vect[1]
    if e_vect[2][1] < 0:    # eighen vector 2 is the direction of brain height, should be aligned with the positive y direction
        e_vect[2]=-e_vect[2]
    rot=numpy.array([e_vect[1], e_vect[2], e_vect[0]]).T 
 
    # Construct empty images of the ROI shape.
    K=D_shape[2]*8/9 # number of slices      
    I=D_shape[0]/4+30 # left/right 
    J=D_shape[1]/4-6 # up/down 
    C_image_data=numpy.zeros((I,J,K),numpy.int16)
    #C_mask_data=numpy.zeros((I,J,K),numpy.uint8)

    # Make sure the the origin voxel index coincide with the origin of World coordinate.
    for k in range(-K/2, K/2):
        logger.debug("working on %d th slide", k+K/2)
        z=k*k_spacing 
        (x_bar,y_bar,z_bar)=rot.dot([0,0,z])+center
        (xd,yd,zd)=distort(x_bar,y_bar,z_bar)

        (id,jd,kd,a)=numpy.dot(dicom_affine_inv,[xd,yd,zd,1])
        if kd < 0: # have not reached the (DISTORTED) image's area
            continue
        elif kd > D_shape[2]: # finished every slice 
            break

        for j in range(-J/2, J/2):
            y=j*j_spacing
            (x_bar,y_bar,z_bar)=rot.dot([0,y,z])+center
            (xd,yd,zd)=distort(x_bar,y_bar,z_bar)
            (id,jd,kd,a)=numpy.dot(dicom_affine_inv,[xd,yd,zd,1])
            if jd < 0:
                continue
            elif jd >= D_shape[1]:
                break

            for i in range(-I/2, I/2):
                # convert the voxel position to the scanner position.
                # currently, just to make the new image have the same origin as the old image, and no tilting
                x=i*i_spacing
                (x_bar,y_bar,z_bar)=rot.dot([x,y,z])+center
                (xd,yd,zd)=distort(x_bar,y_bar,z_bar)
                (id,jd,kd,a)=numpy.dot(dicom_affine_inv,[xd,yd,zd,1])
        
                if id < 0:
                    continue
                elif id > D_shape[0]:
                    break 
                else:
                    C_image_data[i+I/2,j+J/2,k+K/2]=numpy.int16(interpolate_tl(D_image_data,id,jd,kd))
                    assert C_image_data[i+I/2,j+J/2,k+K/2] >= 0
                    #C_mask_data[i+I/2,j+J/2,k+K/2]=numpy.uint8(interpolate_nn(D_mask_data,id,jd,kd))


    # construct nifti image
    hdr=D_hdr.copy()
    hdr.set_data_shape((I,J,K))
    aff=numpy.array([[-i_spacing,0,0,-I/2*i_spacing],[0,-j_spacing,0,-J/2*j_spacing],[0,0,k_spacing,K/2*k_spacing],[0,0,0,1]])
    hdr.set_qform(aff, code='scanner')
    hdr.set_sform(aff, code='scanner')

    img=nibabel.Nifti1Image(C_image_data,affine=aff,header=hdr)
    img.to_filename(C_image_file) 

if __name__=="__main__" :
    import logging
    logger = logging.getLogger("correction")
    logger.setLevel(10)
    if len(logger.handlers) == 0 :
        formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s")
        ch = logging.StreamHandler()
        ch.setLevel(10)
        ch.setFormatter(formatter)
        logger.addHandler(ch)
    if len(sys.argv)==5:
        action=sys.argv[1]
        brain=sys.argv[2]
        domain=sys.argv[3]
        codomain=sys.argv[4]
        globals()[action](brain,domain,codomain)

